/*******************************************************************************

                                DISTANCE

  Title:        Calculate distances among files of electron/pixel.    
  Author:       Hanna Szoke
  Date:         29/12/94
  Function:     This program reads M (in range 2 - 9) sets of binary files 
		representing (unregridded) electron/pixels.  For each 
		pair, it calculates and reports 

		* the rms fractional distances among them;
		* the linear absolute fractional distances among them; and
		* the correlation coefficients among them.
		
                Input parameters are read locally. UTIL.C has worker functions.

*******************************************************************************/

#include "util.h"
#include "dims.h"

static	char	dfilename[MAXNFILES][MAXSTRING] ;
static	real	*dnump[MAXNFILES] ;	/* electron densities */
static	int	M ;			/* number of Np file sets to measure */

static  int	dNhr[MAXNFILES] ;
static	real	*hrNump[MAXNFILES] ;
static	real	*hrNx[MAXNFILES] ;
static	real	*hrNy[MAXNFILES] ;
static	real	*hrNz[MAXNFILES] ;

static	void	do_cc() ;
static	void	do_lin() ;
static	void	do_sparseD() ;
static	void	do_rms() ;
static	void	do_overlap() ;
static	real	rms_fracdist(real *, real *, real *, real *) ;
static	real	lin_fracdist(real *, real *, real *, real *) ;
static	real	sparseD(real *, real *, real *, real *) ;
static	real	corr_coef(real *, real *, real *, real *) ;
static	real	overlap(real *, real *, real *, real *) ;
static	void	printMatrix(real (* )()) ;
static  char	is_sparse(real *) ;


void	distance_main(int argc, char *argv[])
{
     char	list_filename[MAXSTRING] ;
     void	readEpixFile(char *, real *) ;
     int	hrGetLength(char *) ;
     void	hrRenameList(real **, real **, real **) ;
     void	hrReadLists(char *, real *) ;		
     void	hrCompareIndices(real *, real *, real *) ;

     int	m ;			/* counters of files */

  /***********************************************************
   Check for right ballpark, say hello, identify file names.
  ***********************************************************/

     M = argc - optind - 2 ;

     if (M > MAXNFILES) 
	  Fatal("Too many files!") ;
     else if (M >= 2)
          for (m = 0; m < M; m++) 
               sprintf(dfilename[m], "%s", argv[optind+m+2]) ;

     else {
	  prompt("\nPlease enter 2 - 8 file names: ") ;
          sscanf(terminp, "%s %s %s %s %s %s %s %s", dfilename[0], dfilename[1],
	         dfilename[2], dfilename[3], dfilename[4], dfilename[5], 
	         dfilename[6], dfilename[7]) ;

          M = m = 0 ;
          while ((int)strlen(dfilename[m]) > 0)   {
	       M++; 
	       m++;
          }
          if (M < 2) 
	       Fatal("Too few files!") ;

	  sprintf(message, "Running %s on %d files", caller, M) ; 
	  printTwice(message) ;
     }

     hello(caller) ;

  /************************************************************
  Fetch input conditions, set up arrays in physical space.
  Read the binary el/voxel maps. 
  ************************************************************/

     readBasicInput() ;

     for (m = 0; m < M; m++) 
          dnump[m] = (real *) e_malloc(Nptotal*sizeof(real), "distance") ;

     for (m = 0; m < M; m++) {  
          readEpixFile(dfilename[m], dnump[m]) ;

          if (HighRes) {
               sprintf(list_filename, "%s.list", dfilename[m]) ;
               dNhr[m] = hrGetLength(list_filename) ;	
	       if (m > 0) 
	            if (dNhr[m] != dNhr[0])
	                 Fatal(
               "Your high-res list files do NOT have the same length!") ;

               hrNump[m]= (real *)e_malloc(dNhr[m]*sizeof(real), "distance") ;
               hrReadLists(list_filename, hrNump[m]) ;	
	       hrRenameList(&hrNx[m], &hrNy[m], &hrNz[m]) ;

	       if (m > 0) 
                    hrCompareIndices(hrNx[m], hrNy[m], hrNz[m]) ;
          }
     }
  /************************************************************
  Compute and report distances.
  ************************************************************/

     do_rms() ;
     do_lin() ;
     do_sparseD() ;
     do_cc() ;
     do_overlap() ;

}

void	do_rms() 
{
  /************************************************************
  Compute and report rms distances.
  ************************************************************/

     printTwice(
       "\nThe rms fractional distance between 2 el/pix files - A and B - is\n") ;
     printTwice( "\tsqrt[Sum[(a - b)^2]] / (sqrt[Sum[a^2]] + sqrt[Sum[b^2]])\n ") ;
     printTwice( "where a and b are data values from files A and B") ;
     printTwice( "and Sum[] is a sum over all such elements.\n") ;

     if (M == 2) {
          sprintf(message, 
	  "The rms fractional distance between data in %s and %s is %g\n",
          dfilename[0], dfilename[1], 
	  rms_fracdist(dnump[0], dnump[1], hrNump[0], hrNump[1]));
          printTwice(message) ;
     }
     else {
	  printTwice("Rms fractional distances among the files.\n") ;

          printMatrix(rms_fracdist) ;

     }
}
real	rms_fracdist(real *targ1, real *targ2, real *hr1, real *hr2) 
{
     real 	*ptarg1 ;
     real 	*ptarg2 ;
     int 	n ;
     real	sumsqdiff = 0 ;
     real	sumsqa = 0, sumsqb = 0 ;
     real	a, b,  diff ;

     ptarg1 = targ1 ;
     ptarg2 = targ2 ;

     for (n = 0; n < Nptotal; n++) {
          a = *(ptarg1++) ;
          b = *(ptarg2++) ;
	  diff = a - b ;
	  sumsqdiff += diff * diff ;
	  sumsqa    += a * a ;
	  sumsqb    += b * b ;
     } 
     if (HighRes) {
	  for (n = 0; n < dNhr[0]; n++) {
               a = *(hr1++) ;
               b = *(hr2++) ;
	       diff = a - b ;
	       sumsqdiff += diff * diff ;
	       sumsqa    += a * a ;
	       sumsqb    += b * b ;
          }
     }
     return (sqrt(sumsqdiff) / (sqrt(sumsqa) + sqrt(sumsqb))) ;

}

void	do_lin() 
{
  /************************************************************
  Compute and report linear absolute fractional distances.
  ************************************************************/

     printTwice( 
  "\nThe linear absolute fractional distance between 2 el/pix files - A and B - is\n") ;
     printTwice( "\tSum[|a - b|] /  Sum[a + b]\n ") ;
     printTwice( "where a and b are data values from files A and B") ;
     printTwice( "and Sum[] is a sum over all such elements.\n") ;

     if (M == 2) {
          sprintf(message, 
	  "The linear abs. fractional distance between data in %s and %s is %g\n",
          dfilename[0], dfilename[1], 
	  lin_fracdist(dnump[0], dnump[1], hrNump[0], hrNump[1])) ;
          printTwice(message) ;
     }
     else {
	  printTwice(
        "The linear absolute fractional distances among the files.\n") ;

          printMatrix(lin_fracdist) ;
     }
}
real	lin_fracdist(real *targ1, real *targ2, real *hr1, real *hr2) 
{
     real 	*ptarg1 ;
     real 	*ptarg2 ;
     int 	n ;
     real	sum = 0, diff = 0 ;

     ptarg1 = targ1 ;
     ptarg2 = targ2 ;

     for (n = 0; n < Nptotal; n++) {
	  diff += fabs(*(ptarg1) - *(ptarg2)) ;
	  sum += (*(ptarg1++) + *(ptarg2++)) ;
     } 
     if (HighRes) {
	  for (n = 0; n < dNhr[0]; n++) {
	       diff += fabs(*(hr1) - *(hr2)) ;
	       sum += (*(hr1++) + *(hr2++)) ;
          }
     }
     return (diff / sum) ; 

}

void	do_sparseD() 
{
  /** Compute and report sparse distances. **/

     printTwice( "\nThe distance between 2 sparse el/pix files - A and B - is\n") ;
     printTwice( "\tSum[|a - b|^2 / (a + b)] / Sum[(a + b)]]\n ") ;
     printTwice( "where a and b are data values from files A and B") ;
     printTwice( "and Sum[] is a sum over all such elements") ;
     printTwice( "for which (a + b) is > 1% of (a_max+b_max).\n") ;
     printTwice( "It is assumed that the signals have Poisson random noise.") ;
     printTwice( "The measure should not be dominated by small values of the density.") ;

     if (M == 2) {
          sprintf(message, 
	  "The distance between data in %s and %s is %g\n",
          dfilename[0], dfilename[1], 
	  sparseD(dnump[0], dnump[1], hrNump[0], hrNump[1])) ;
          printTwice(message) ;
     }
     else {
	  printTwice("The distance among the files.\n") ;
          printMatrix(sparseD) ;
     }
}
real	sparseD(real *targ1, real *targ2, real *hr1, real *hr2) 
{
     real 	*ptarg1, *ptarg2 ;
     int 	n ;
     real	sum, diff ;
     real	sumsqdiff = 0,	sumsum = 0 ;
     real	ar_max1 = -10000., ar_max2 = -10000. ;
     real	epsilon ;

     for (n = 0, ptarg1=targ1, ptarg2 = targ2; n < Nptotal; n++, ptarg1++, ptarg2++) {
	  if (*(ptarg1) > ar_max1)
	       ar_max1 = *(ptarg1) ;
	  if (*(ptarg2) > ar_max2)
	       ar_max2 = *(ptarg2) ;
     }
     epsilon = (ar_max1 + ar_max2) / 100. ;

     for (n = 0, ptarg1=targ1, ptarg2 = targ2; n < Nptotal; n++, ptarg1++, ptarg2++) {
	  diff = fabs(*(ptarg1) - *(ptarg2)) ;
	  sum = (*(ptarg1) + *(ptarg2)) ;
	  if (sum > epsilon) {
	       sumsqdiff += diff * diff / sum ;
	       sumsum += sum ;
          }
     } 
     if (HighRes) {
	  for (n = 0; n < dNhr[0]; n++) {
	       diff = fabs(*(hr1) - *(hr2)) ;
	       sum = (*(hr1++) + *(hr2++)) ;
	       if (sum > epsilon) {
	            sumsqdiff += diff * diff / sum ;
	            sumsum += sum ;
               }
          }
     }
     return (sumsqdiff / sumsum) ; 
}

void	do_overlap() 
{
     int	m ;
     char	spflag ;

     /* Are files sparse enough to warrant this measure? */

     for (m = 0; m < M; m++) {
	  spflag = is_sparse(dnump[m]) ;
	  if (!spflag) {
	       printTwice(
	       "Not all the data sets are sparse; overlap report is skipped") ;
	       return ;
          }
     }

  /** Compute and report measure of overlap. **/

     printTwice( "\nThe overlap between 2 el/vox files measures the extent") ;
     printTwice( "to which comparable entries have significant electron counts,") ;
     printTwice( "but doesn't compare the values of these counts.") ;

     if (M == 2) {
          sprintf(message, "\nThe overlap between data in %s and %s is %g\n",
          dfilename[0], dfilename[1], 
	  overlap(dnump[0], dnump[1], hrNump[0], hrNump[1])) ;
          printTwice(message) ;
     }
     else {
          printTwice(
	  "Note: matrix is not symmetric! 'ave' measures ROW average.\n") ;
	  printTwice("The overlap among the files follows.\n") ;
          printMatrix(overlap) ;
     }
}
char	is_sparse(real *targ)
{
     real 	*ptarg ;
     int 	n, count ;
     real	ar_max = -10000. ;


     for (n = 0, ptarg=targ; n < Nptotal; n++, ptarg++) 
	  if (*(ptarg) > ar_max)
	       ar_max = *(ptarg) ;
     
     for (n = 0, ptarg=targ, count = 0; n < Nptotal; n++, ptarg++) 
	  if (*(ptarg) > 0.01 * ar_max)
	       count++ ;
     
     return (count < Nptotal/100) ;
}
real	overlap(real *targ1, real *targ2, real *hr1, real *hr2) 
{
     real 	*ptarg1, *ptarg2 ;
     int 	n, count1, count2 ;
     real	ar_max1 = -10000., ar_max2 = -10000. ;

     for (n = 0, ptarg1=targ1, ptarg2=targ2; n < Nptotal; n++, ptarg1++, ptarg2++) {
	  if (*(ptarg1) > ar_max1)
	       ar_max1 = *(ptarg1) ;
	  if (*(ptarg2) > ar_max2)
	       ar_max2 = *(ptarg2) ;
     }
     for (n = 0, ptarg1=targ1, ptarg2=targ2, count1 = count2 = 0; 
	  n < Nptotal; n++, ptarg1++, ptarg2++) 

	  if (*(ptarg2) > 0.005 * ar_max2) {
	       count2++ ;
     
	       if (*(ptarg1) > 0.005 * ar_max1)
	           count1++ ;
          }
	  
     return (count1/(real)count2) ; 
}

void	do_cc() 
{
  /************************************************************
  Compute and report correlation coefficients.
  ************************************************************/

     printTwice(
     "\nThe correlation coefficient between 2 el/pix files - A and B - is\n") ;
     printTwice( "\tSum[del(a) * del(b)] / sqrt[Sum[del(a)^2] * Sum[del(b)^2)]]\n") ;

     printTwice( "where del(a) is a - a_ave and del(b) is b - b_ave") ;
     printTwice( "and Sum[] is a sum over all such elements.\n") ;

     if (M == 2) {
          sprintf(message, 
	  "The correlation coefficient between data in %s and %s is %g\n",
          dfilename[0], dfilename[1], 
	  corr_coef(dnump[0], dnump[1], hrNump[0], hrNump[1])) ;
          printTwice(message) ;
     }
     else {
	  printTwice( "The correlation coefficients among the files.\n") ;
          printMatrix(corr_coef) ;
     }
     return ;
}
real	corr_coef(real *targ1, real *targ2, real *hr1, real *hr2) 
{
     real 	*ptarg1 ;
     real 	*ptarg2 ;
     real	mean1 = 0 ;
     real	mean2 = 0 ;
     real	num = 0 ;
     real	denom1 = 0, denom2 = 0 ;
     real	d1, d2 ;
     int 	n ;
     int	N ;

     for (n = 0, ptarg1 = targ1; n < Nptotal; n++, ptarg1++) 
	  mean1 += *ptarg1 ;

     N = Nptotal; 

     if (HighRes) {
	  for (n = 0; n < dNhr[0]; n++) 
	       mean1 += *(hr1+n) ;
          N += dNhr[0] ;
     }
     mean1  *= 1./(float) N ;

     for (n = 0, ptarg2 = targ2; n < Nptotal; n++, ptarg2++) 
	  mean2 += *ptarg2 ;

     if (HighRes) 
	  for (n = 0; n < dNhr[0]; n++) 
	       mean2 += *(hr2+n) ;

     mean2  *= 1./(float) N;

     for (n = 0, ptarg1 = targ1, ptarg2 = targ2; n < Nptotal; 
	  n++, ptarg1++, ptarg2++) {

	  d1 = *ptarg1 - mean1 ;
	  d2 = *ptarg2 - mean2 ;
	  num += d1 * d2 ;
	  denom1 += d1 * d1;
	  denom2 += d2 * d2 ;
     }
     if (HighRes) 
          for (n = 0; n < dNhr[0]; n++) {

	       d1 = *(hr1+n) - mean1 ;
	       d2 = *(hr2+n) - mean2 ;
	       num += d1 * d2 ;
	       denom1 += d1 * d1;
	       denom2 += d2 * d2 ;
          }

     return (num / sqrt(denom1*denom2) ) ;
}

void printMatrix(real (* distfunc)(real *, real *, real *, real *))
{
     int	m, m1, m2 ;		/* counters of files */
     int	l, len ;		/* counters of length message buffer */

     real	average = 0 ;
     real	maver[MAXNFILES] ;
     real	value;			/* value returned by rms_fracdist(() */
     char	alpha[MAXNFILES] ;

  /************************************************************
  Compute and report distances.
  ************************************************************/

	  for (m = 0; m < M; m++) 
	       alpha[m] = 'a' + m ;
	  for (m = 0, l = 0, len = 0; m < M; m++, l+=len) 
	       len = sprintf(message+l, "\t%5s%c", "     ", alpha[m]) ;
          printTwice(message) ;

	  for (m1 = 0; m1 < M; m1++) {
	       len = sprintf(message, "%5s%c", "     ", alpha[m1]) ;
	       for (m2 = 0,l=len, maver[m1] = 0; m2 < M; m2++,l+=len) {
		    value = distfunc
			 (dnump[m1], dnump[m2], hrNump[m1], hrNump[m2]) ;
		    if (m2 > m1) 
			 average += value ;
		    if (m2 != m1) 
			 maver[m1] += value ;
		    len = sprintf(message+l, "\t%6.4f", value) ;
               }
               printTwice(message) ;
          }

          len = sprintf(message, "\n   ave") ;

       	  for (m2 = 0,l=len; m2 < M; m2++,l+=len) {
	       len = sprintf(message+l, "\t%6.4f", maver[m2] / (M-1)) ;
          }
          printTwice(message) ;

	  sprintf(message, "\nGrand average = %g\n", average / (M*(M-1)/2)) ;
          printTwice(message) ;

	  printTwice("Where:") ;
	  for (m = 0; m < M; m++) {
	       sprintf(message, "\t%c%s%s", 
		    alpha[m], " stands for ", dfilename[m]) ;
               printTwice(message) ;
          }
}
